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Abstract 

We introduce a family of symplectic, linearly-implicit and stable integrators for 
mechanical systems. When used in conjunction with penalty methods (i.e., methods 
that consist in replacing constraints by stiff potentials), these integrators acceler- 
ate the numerical simulation of mechanical systems with holonomic constraints by 
employing coarse timesteps and bypassing the resolution of nonlinear systems. Al- 
though penalty methods are well known and widely employed, a general and rigorous 
proof of their accuracy appeared to be lacking; such a proof is also provided in this 
paper. 

1 Introduction and main results 

In this paper, we introduce a new family of symplectic, linearly-implicit and stable 
integrators for the following mechanical dynamics 

(1) 

-VVix) ^ ' 

where V G C^(M'^) is a potential energy function and M & n x n constant, symmetric 
mass matrix. This new family of integrators (SyLiPN) is parameterized by a constant 
j3 and is of the form 

Integrator 1. Symplectic Linearized Push- forward Newmark (SyLiPN): 

^ Xk+i = Xk + hvk + \h^ak 
Vk+i =Vk + \h{ak + flfc+i) (2) 
au = -M^^W{xk) - M-iHessF(xfc)/3/i2afc 



Remark 1.1. Notice that the third line, i.e., the force evaluation, can be rewritten as 

Ofc = -(M + Hessy(xfe)/3/i2)-ivy(xfc) (3) 
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This evaluation, however, should be executed by solving a symmetric linear system 
instead of inverting a matrix for computational efficiency. In this sense, SyLiPN is 
linearly implicit. 

Remark 1.2. We refer to Integrator [J] for a simplification of SyLiPN when V can be 
decomposed as the sum of a soft and a stiff potential. 

In Section [21 we will show how SyLiPN is obtained by linearizing the Push-forward 
of the Newmark family of implicit integrators |24^ l30l |22] . We will also prove (in Section 
[2]) that SyLiPN has the following properties. 

Theorem 1.1. SyLiPN (Integrator\^ : 

1. Is linearly unconditionally stable if f3 > 1/4. 

2. Has a 3rd order local error. 

3. Is symplectic. 

Our motivation for developing SyLiPN originates from (but is not limited to) a need 
for such integrators in penalty methods, i.e., methods based on replacing holonomic 
constraints with stiff potentials [331 EH [26] as described below (we note that SyLiPN 
is also relevant to applications, such as molecular dynamics, where stiff potentials are 
replaced with holonomic constraints [71 [28] ) . 

Let ^ 

S{q{t)) := ^ lq{tfMq{t) - V{q{t)) dt (4) 

be the action functional associated with the mechanical dynamical system ([!]). The 
evolution of this mechanical system under the holonomic constraint g{q) = can be 
obtained as the critical trajectory on the constrained manifold, i.e., the solution of: 

5S/5q = and for ah t, q{t) € c/~^(0) (5) 

This constrained dynamic can also be shown as the solution to the following Differential 
Algebraic Equation system: 

{Mq = -VV{q) + \^Vg{q) 
\9{q) =0 

The idea behind penalty methods is to replace the rigid constraints by stiff poten- 
tials that enforce the holonomic constraints up to small deviations. More precisely, the 
potential energy V{q) is modified to be V{q) + ^oj'^giq)^ g{q), and the solution of the 
constrained system ([5|) is approximated by the solution of the following unconstrained 
mechanical system. 

Mq = -VV{q) - oo^g{qfVg{q) (7) 

where lo is chosen to be large enough for accuracy. Paraphrasing [26], the problem with 
this approach is that "as a result of stiffness, the numerical differential equation solver 
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takes very small time steps, using a large amount of computing time without getting much 
done^\ With SyLiPN, this issue can now be solved by employing coarse timesteps and 
bypassing the resolution of nonlinear systems. Up to the authors' knowledge, SyLiPN is 
the first symplectic, stable integrator for such systems that does not require the resolution 
of a nonlinear system at each timestep. 

Although penalty methods are well known and widely employed, a rigorous justifica- 
tion of this accuracy appeared to be lacking. Our analysis shows that the solution of ([7|) 
converges to the solution of ([6]) as w — )■ oo in the sense of two-scale flow convergence (see 
Definition 1.1 in [31j). More precisely, we will prove the following theorem in Section 
[3] (throughout this analysis, 'bounded' means having a norm no larger than a constant 
independent of w). 

Theorem 1.2. Denote by q^{t) the solution to ^ with g'^(O) = and g^(0) = go 
(where g{qo) = and ^g{qo) = V(7(go) ' Qo = 0)- Suppose that M is non-singular, V{-) 
is bounded from below, V{q) diverges towards infinity as \q\ — )• oo, V{-) and g{-) are 
with hounded derivatives, and for all q G g~^{Q), Vg{q) has a constant rank equal to the 
codimension of the constrained manifold, then 

X{t) := - hm lim - / u^g{q^{s))ds (8) 

exists. Also, the solution to 

'um = -VV(q{t)) + MtfVg(g{t)) 

g(0) = qo 

m = qo 

which we denote by q{t) (q{t) is also the solution to exists and satisfies 

q^ 4 q (10) 

in the sense that for all bounded t > and all hounded and uniformly Lipschitz continuous 
test function ip, 

lim lim - / ^{q^{s)) - ^{q{s)) ds = (11) 



In sectional we will investigate three numerical examples: (1) double pendulum (2) a 
chain of many pendulums, i.e., an approximation of continuous ropes, hairs and soft tree 
branches (these examples share the same constraints, differ in the soft potentials, and 
are of significant importance in computer graphics) (3) the dynamics of water cluster, 
i.e., multiple water molecules that interact with each other. Numerical results show the 
accuracy of SyLiPN and significantly improved computation time for all three examples. 
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1.1 On implicit methods, penalty methods and constrained dynamics 



Implicit methods are widely used for stable integrations, and the symplectic members of 
them have been successful in long time integrations of mechanical systems [9j. Implicit 
methods require solving nonlinear systems, which oftentimes cost more than one iter- 
ation in a nonlinear solver (e.g., Newton's method [25]) when symplecticity is desired. 
Indeed, for symplectic integrators, if one solves the nonlinear updating equations par- 
tially (for instance, by carrying out only the first step in a Newton's solver) or linearizes 
the nonlinearity (which can be shown to be equivalent to the former), in most cases sym- 
plecticity is lost. Such a loss results in unsatisfactory long time performances, including 
drifts in energy and momentum (imagine a tennis video game in which the ball does 
not bounce back from the racquet constrained in the player's grip), or even instability. 
SyLiPN is linearly-implicit in the sense that it only requires the resolution of a linear 
system at each integration step, which can be computed efficiently by methods provided 
in, for instance, [2] and references therein. SyLiPN inherits both the stability and the 
symplecticity from Newmark, however at a much lower computational cost. In other 
words, the choice of Push-forward Newmark in SyLiPN allows the nonlinear implicit 
updates to be solved by the first Newton iteration without losing symplecticity. 

As an important application of SyLiPN, mechanical systems with holonomic con- 
straints can be simulated more efficiently than by existing methods, which include 
generalized coordinate approaches (e.g., p^J) and Lagrange multiplier methods (e.g., 
[271 E Uni [Ml |T7|) — both are implicit in general cases. For doing so, we first approxi- 
mate the constrained dynamics (|5|) by the modified system (|7|) in which rigid algebraic 
constraints are replaced by stiff nonlinear springs, and then employ SyLiPN to integrate 
the modified system ([7]) with a large timestep (which does not require the resolution of 
the stiffness in the equations). The strategy that consists in replacing holonomic con- 
straints by stiff potentials (and vice versa) is, by now, widely acknowledged and used. In 
molecular dynamics, for instance, stiff oscillatory molecular bonds are replaced with rigid 
bond length constraints for the purpose of allowing large timesteps (e.g., fTiiSSj). The 
reverse point of view, which corresponds to replacing rigid bond length constraints with 
stiff oscillatory bonds, is of significant and practical importance in computer graphics 
and related fields (we refer for instance to [Ml (Ml [20] ) . 

Traditional methods for simulating the constrained mechanical system defined by ([5|) 
include (but are not limited to): introducing generalized coordinates on the constrained 
manifold (e.g., [13j), or using Lagrange multipliers (e.g., SHAKE [27], RATTLE [T], SET- 
TLE [LQI, LINCS [23], M-SHAKE [17]). The equivalence between these two approaches 
is now well-established (see for instance [35] for a discussion of the mathematical justi- 
fication). These numerical methods (e.g., [131 [271 [D HOl ISS HZ]) allow a o(l) integration 
step, but they also require solving nonlinear systems at each step (which significantly 
slows down computations). Another strategy would be to use a classical symplectic in- 
tegrator (such as Velocity Verlet) to explicitly simulate the modified mechanical system 
([7|) at the cost of using o{u}^^) timesteps. The best strategy depends on the system: if 
the constraints are high dimensional, an explicit integration, even with small steps, will 
still be faster, because the nonlinear solve is too expensive; however, if there are just a 
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few constraints, generalized coordinate and Lagrange multiplier methods will be more 
efficient. 

SyLiPN (Integrator [1]) allows for large (o(l)) integration timesteps with the limited 
cost of a linear solve per iteration (instead of a nonlinear one). SyLiPN also remains 
accurate with o(l) timestep when applied to ([7]) because the stiffness in this system 
results in a fast dynamics that (locally) converges to a point distribution (see analysis in 
Section [3]), whose effective contribution to the slow dynamics can therefore be captured 
by an implicit method [21j . Therefore, SyLiPN will be both efficient and accurate for 
constrained dynamics. Let us also observe that, although efficient, M-SHAKE is 
limited to systems with distance constraints, and SyLiPN does not suffer from this 
restriction: it is a more general approach that can be applied to arbitrary holonomic 
constraints. 



2 SyLiPN: a family of symplectic, linearly-implicit and sta- 
ble integrators 

The following Newmark family of algorithms have been widely used in structure dynam- 
ics [M]: 

Integrator 2. Newmark: 

' Qk+i = Qk + Hk + ^[{l - 2l3)ak + 2Pak+i\ 

< qk+i =qk + h[{l - 'y)ak + TOfe+i] (12) 
ak =M-\-VV{qk)) 



It was known [22] that Newmark is 2nd-order accurate when 7 = 1/2 and 1st- 
order otherwise, and it is generally implicit when /3 7^ 0. It was also shown [15] that 
Newmark is variational for arbitrary /3 when 7 = 1/2 (we will restrict ourselves to this 
case throughout this paper). However, it is worth noticing that the symplectic form that 
Integrator [2] preserves is not the canonical one. In fact, it was shown \6^\ I22| that if one 
pushes forward the Newmark integrator by a coordinate transformation 77 : TQ — )• TQ 
defined as 

{x,v):=r^{q,q) = {q + Ph^M-^VV{q),q) , (13) 
then we obtain an integrator that preserves the canonical symplectic form on T*Q: 

Integrator 3. Push-forward Newmark: 

Xk+i = Xk + hvk + \h'^ak 

Vk+i =Vk + \h{ak + ak+i) (14) 
ak =-M-^VV{xk + Ph^ak) 
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The Newmark and Push-forward Newmark schemes can be shown to be uncondition- 
ahy hnearly stable if /? > 1/4 [5l[30], and Newmark was further shown to be nonhnearly 
stable in the same case under several assumptions [12]. In addition, Push-forward New- 
mark was shown to be stable near stable fixed points in general nonlinear settings, unless 
specific resonances occur |29) . Nevertheless, there are nonlinear cases in which Newmark 
is no longer stable [6l[19]. In fact, few convergent methods are unconditionally stable 
for arbitrary nonlinear systems to the authors' knowledge (e.g., see a discussion in [37J ) . 

Now, linearize the nonlinear force in Push-forward Newmark at each step by Taylor 
expansion at x^, so that a nonlinear implicit equation becomes a linear one. One obtains 
Integrator [U i.e., SyLiPN. 

Remark 2.1. Since Newmark or Push- forward Newmark requires solving a nonlinear 
system at each step, SyLiPN exhibits a speed advantage. We will quantify this advantage 
numerically in Section HI 

Next, we will show that SyLiPN is stable, convergent, and symplectic. SyLiPN is, 
so far, the only method with these properties that does not require the resolution of a 
nonlinear system at each coarse timestep. 

Theorem 2.1 (Stability). SyLiPN (Integrator [1^ is linearly unconditionally stable if 
P > 1/4. 

Proof. The definition of linear unconditional stability for mechanical system is that when 
V{x) = ^x^Kx with K positive definite, the integrator produces bounded trajectories 
for all time for all integration step h. 

Plotting V into the SyLiPN updating rule ([2]), we obtain 

Xk+i = Xk + hvk + ^/i^Ofc 
' Vk+i = Vk + iH^^k + flfc+i) (15) 
^Mok = —Kxk — Kj3h?ak 

Notice that, since M and K are positive definite, there exists an invertible matrix Q, 
such that Q^MQ = I and Q^KQ = A for a diagonal matrix A (this is known as 
simultaneous diagonalization by congruence; see Corollary 1.7.18 in [11]). Therefore, 
letting Oi = Q~^ai, xi = Q~^Xi, and Vi = Q~^Vi, we have 

Xk+i = Xk + hvk + ^/i^flfc 

Vk+i =Vk + \h{ak + dk+i) (16) 
MQdk = -KQxk - KQl3h^dk 

Multiplying the last equation by on the left, we get 

dk = -Axk - Aph^dk (17) 

Therefore, one can assume without loss of generality that M is the identity and K is 
a diagonal matrix; since the following analysis applies to diagonal elements, we can, 
without loss of generality, further assume M and K to be scalars. 
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Now, from the 3rd equation in (|15|) . we have 

ak = -{M + K^h^y^KXk (18) 

Substituting the above expression for at and Uk+i into the updating rule, we obtain 

I-ah'^/2 
ah{ah'^/A-I) I 



Xk+l 









hi 





Xk 




Vk_ 



(19) 



where a = (M -\- K f3h'^) ^K. Algebraic computations show that the updating matrix in 
has the following eigenvalues: 



A 



1,2 



(2 - ah^ ± hyja'^h'^ - 4a) 



(20) 



It is not difficult to show that if a^h?' — 4a < then ||Ai^2|| = 1 (and hence the trajectory 
remains bounded). Since M > and X > 0, we only need to show that ah? — 4 < 0. 
Because 

Kh'' , Kh'' 



ah 



< 



M + K/3h'^ - Kj3h'^ 
if /3 > 1/4, the above is < 4. This finishes the proof. 
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(21) 
□ 



Remark 2.2. This theorem can be easily understood: since only linear test problems 
are considered, SyLiPN is identical to Push-forward Newmark, which is equivalent to 
Newmark in terms of stability. The stability statement follows because Newmark is 
well-known to be unconditionally stable when j3 > 1/4. 

Remark 2.3. For the special case of a stiff potential V{x) = Vo(x) + e~^V\{x), the 
following modification of SyLiPN will also be unconditionally linearly stable as long as 
/? > l/4 + 0(e): 

Integrator 4. Simplified SyLiPN for stiff systems (e » 1): 



Xk+i = Xk + hvk + ^h'^ttk 
Vk+i =Vk + \h{ak + ak+i) 

ak = -M-i(Vyo(^fc) + e-^VVi{xk)) - M-^e-^YiQssVi{xk)l5K^ 



(22) 



a.k 



i.e., the Hessian of the small Vq can be neglected. This is because assuming Vq{x) 
^x^Kqx and Vi{x) = |: 



\x^Kox and Vi{x) = ^x^Kix, then a in (fT9]) will be 



" " M + e-iKi/3/i2 

It can be shown that the stability condition a^/i^ — 4a < is satisfied for all h as long 
as P > 1/4 + 0(6), because ([23]) can be rewritten as 

Ko/i^ - 4M < e"^i^i(4/3 - 1)/!^ (24) 
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and if this were to violate, the only possibility is that the right hand side is small enough 
because h = 0{^/e), but then the left hand side is < 0, and hence this hypothetical 
violation will not happen. 

Ase— )-0, /3>l/4 will be enough to guarantee stability. 

Remark 2.4. [12] uses an energy bound based on an assumption to demonstrate the 
nonlinear stability of Newmark. The same energy bound applies to SyLiPN, because 
Newmark and SyLiPN differ in force estimations by only a high order term of 0{h'^), 
which will not affect the leading term in the energy bound. The assumption introduced 
there, however, can not be checked a priori for either Newmark or SyLiPN. Possible 
violations of this assumption may result in a nonlinear instability of Newmark or SyLiPN. 
As commented before, unconditional stability for arbitrary nonlinear systems is beyond 
the scope of current researches. 

Remark 2.5. For possible improvements on nonlinear stability, one may resort to lin- 
earizations of more stable methods (such as those in il8j). However, few methods remain 
symplectic after the linearization. 

Theorem 2.2 (Consistency). Consider an integrator for Q given by the following 
updating rule: 

Xk+i = Xk + hvk + \h^ak 

Vk+i =Vk + \h{ak + Ofe+i) , (25) 

^ ak = -M-i Vy (xfc) - M-^f{xk)h^ak 

where f S C{Q) is an arbitrary function. 

IfV^ C'^{Q), this integrator has a 3rd order local error, i.e., start with Xk,Vk at time 
kh, and denote by Xk+i,Vk+i the exact solution at time {k + l)h and by Xk+i,Vk+i the 
numerical solution after one-step update, then Xk+i — Xk+i = 0{h'^) and Vk+i — ffc+i = 
0{h^). 

Proof. Assume without loss of generality that M = I. Writing a(-) = — Vy(-), we have 
ak = a(xfc)/(l + fixk)h^) = a(xfc) + Oih^) (26) 

and 

a{xk+i) = a{xk)-^{xk+i-Xk)a' {xk)-^0 {{xk+i - Xkf) = a{xk)+hvka' {xk)+0{h^) (27) 
Since the exact dynamics ([TJ is governed by 

(28) 
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if smoothness of the solution is assumed, we obtain by Taylor expansion that 



Xk+i = Xk + hvk + yUfc + 0{h ) 

= xk + hvk + —ak + OQv^) = Xk+i + 0{h?) 
Vk+1 = Vk + ha{xk) + —a{xk)vk + 0{h^) 

= Vk + ^{a{xk) + a{xk) + ha'{xk)vk) + 0{h^) 

= Vk+'^iak + Oik") + a{xk+i) + 0{h^)) + 0{h^) 

= Vk + ^{ak + Oik") + ak+i + Oik")) + 0{h^) = Vk+i + Oih^) 



(29) 



(30) 
□ 



Corollary 2.1. SyLiPN (Integrator\^ has a 3rd order local error. 

Proof. This is because SyLiPN (l2|) is a special case of (f25|) with f{xk) = Hessy(xfc)/3. □ 

Remark 2.6. By the celebrated Lax-Richtmyer equivalence theorem [20], a ©(/i^^^) 
local error (also known as 'consistency' when p > 1) together with stability will lead to 
a 0{W) global error (i.e., 'convergence' for p > 1). For our case, if SyLiPN is stable, 
then it is 2nd order convergent (with a ©(/i^) global error). 



Theorem 2.3 (Symplecticity). Consider an integrator given by the following updating 
rule: 

Xk+l 

(31) 



Xk + hvk + hh'^ak 



Vk+i =Vk + \h{ak + Ok+i) , 
[^ttk = fi^k) 

where f £ C^{Q) is an arbitrary function. This integrator is symplectic. 
Proof. Compute the Jacobian of Eq. ([3T]) (viewed as Xk,Vk ^ Xk+i,Vk+i): 



V :-- 



(9a; J; 
dvk+i 



dvi. 



I+^h'f'{xk) hi 

\h {f(xk) + (/ + ^h^ixk)) f'{xk+i)) I + \h^f'{xk+i) 



(32) 
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Let J 



/ 
-/ 



be the canonical symplectic matrix. Basic algebra shows that 



V^JV 



A 
B C 



where 



'A={I + (/ + ^h^fixk+i)) - ifixk) + + Ih^fixu)) f'ixk+i)) = 

B={I + ^h^fixk)) {-I - ^h^fixk+i)) + i/i^ if'ixk) + (/ + ^h^fixk)) f'{xk+i)) 
C = h{-I- Ih^f'ixk+i)) +h{l+ lh^f'(xk+i)) = 

which is indeed the definition of symplecticity JT) = J . □ 

Corollary 2.2. SyLiPN ( Integrator is symplectic. 

Proof. In SyLiPN, is given by ([3]), which is a function of Xk. Therefore, SyLiPN is in 

the form of (j3ip and hence symplectic. □ 



Remark 2.7. Linearizing Newmark or implicit midpoint will not result in a symplectic 
method, although the resulting integrator will be linearly- imp licit, stable and 2nd-order. 
It is rare that the linearization of an implicit method is symplectic. See also Section [4. II 
for a numerical demonstration on the consequence of loss of symplecticity. 

On the other hand, Theorem 12.31 provides a family of symplectic integrators, but few 
choices of /(•) will result in convergent and stable ones. 

SyLiPN lies in the intersection of linearized approaches and the family of symplectic 
integrators given by ([3T]) . 



Remark 2.8. By Theorems 12.21 and 12. 3|, simplified SyLiPN for stiff systems (Integrator 
HI) is also 2nd-order convergent (if stable) and symplectic. It is also unconditionally 
linearly stable according to Remark [ 



3 Proof of Theorem 11.2 

Consider the solution of d?]). Due to energy conservation in the modified system, 
qMq/2 + V{q) + ^uP' g{qY^ g{q) remains constant in time, and therefore ^uP' g{q)'^ g{q) 
remains bounded independently from oj as long as V{-) is bounded from below (a re- 
quirement that needs to be satisfied only in the domain containing the solution). As a 
consequence, g{q) is at most 0{1/ljj) and the constraint is satisfied approximately. In 
fact, g{q) can be further shown to be of the order of 0(l/a;^) in most cases (we refer to, 
for instance, [16], or Remark 13. II and the multiscale analysis below). 

Theorem 11.21 shows that the modified system ([7]) approximates the Lagrange multi- 
plier approach represented by the DAE ^ and describe the associated topology. Observe 
that there are more than 2000 papers referring to the modified system, i.e., the penalty 
approach [33l |36l |26] . 

The proof of Theorem 11.21 is based on two observations: g{q^) is close to 0, and the 
two scale fiow convergence (F-convergence) holds when the constrained manifold is flat. 
We first elaborate these facts by the following lemmas. 
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Lemma 3.1. // the potential V{-) is bounded from below, then there exists a constant 
C, such that Wgiq'^is))]] < C/uj for all s. Moreover, ifV{q) diverges towards infinity 
as \q\ — )• oo, then the solution is bounded, i.e., there exists a constant C such that 
\\q^{s)\\<C. 

Proof. Notice that the energy [cf^'Mci^ /2 + V{q'^) + uj'^giq'^)'^ g{q'^) in the modified 
system ([7]) is conserved and bounded (due to the initial condition). Therefore, V{-) 
being bounded from below and g^g > imply that oj"^ g{q^)'^ g{q^) = Hence 
g{q^{s)) = 0{l/u:). 

By a similar energy argument, since [(7'^]^Mg'^/2 > and g{q'^)^ g{q^) > 0, V{q'^) is 
bounded from above too, which implies that q^ remains bounded. □ 

Lemma 3.2. Consider the solution to a conserved mechanical system 

= /2(x'^,y")-a;Vy"fV5(r) 

where and y^ are both vectors, and x^iQ) = xq, x^{{)) = xq, ?/'^(0) = yo, y'^{0) = yo- 
Suppose that fi, f2 and Vg are with bounded derivatives, x^ and y^ are bounded, 
5(2/0) = and ^5(2/0) = 0, and g{-) has a non- degenerate Jacobian in a neighborhood of 
yo, then the limit 

A(0:=-lim lim i / ^ u^g{y^{s))ds (34) 
exists and is finite. Denote by x{t), y{t) the solution to 

\x = fi{x,y) 



[y = f2{x,y) + X^Vg{y) 
with the same initial conditions xq, xo,yo,yo, then as w —)• 00, 



(35) 



A y (36) 

U(r) ^ 



Proof. We employ the multiscale averaging theory studied in |31) to demonstrate the 
convergence. Here the constrained manifold is a submanifold with any x and y satis- 
fying g{y) = 0. x'^ is a slow variable and its evolution corresponds to the constrained 
dynamics, y"^ is a fast variable corresponding to a fluctuating deviation from the con- 
strained manifold at a characteristic timescale of 0(1 /uj) , and it lies in the normal bundle 
of the constrained manifold. 

First, consider the case in which g{-) is a linear function g{y) = Cy'^ (the affine case 
can be similarly treated by shifting y). The dynamics of y is governed by 

r = /2(x-,y")-^Vc^C (37) 
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This is a forced harmonic oscillator, and its solution can be written as: 

= r /2(x"(s), s\n{ujCs)C-^/uj ds (38) 

Jo 

where C = \/ C^C is the well-defined matrix square root, and matrix sin is defined either 
by Taylor expansion or diagonalization. Notice that there is no propagation of initial 
conditions, because the initial conditions satisfy the constraint and are therefore null. 

It can be shown from (|38|) (for instance, by Lemma 3.8 in [52]; the idea is that an 
addition 1/uj comes from the sin due to integration by parts) that y^{t) is 0{uj~'^) at 
least up to t = o(l), and y'^ is asymptotically periodic (because ([57|) is asymptotically 
linear) and hence locally ergodic (as uj — t- oo). 

Since is locally ergodic, well-defines A, and Theorem 1.2 in [3T] guarantees 
that the effective equation for (f37|l is 

y = h{x'^,y) + X^C, (39) 

F 

in the sense that y'^ — > y and x'^ — >• x. Notice that the convergence on x is in the strong 
sense, i.e., lim(^_j.oo x^(t) — )• x{t) for all bounded t > 0. This is because x is purely slow, 
for which case, F-convergence implies strong convergence. 

Now consider a fully nonlinear g{-) with a non-degenerate Jacobian. Lemma 13.11 
gives that g{y^) = ©(l/w). Since g{-) has a non-degenerate Jacobian, and y^ is in a 
bounded and therefore compact set, we also have y^ — yo = 0{l/uj). Consequently, the 
dynamics of y^ actually approaches that of a forced oscillator (with equilibrium at yo) 
at a 0(1/0;) timescale, because g{-) approaches a linear function (which is its first order 
Taylor expansion): 

uj''g{yn^Vg{yn 

oo\Vg{yo){y'' - yof + 0{u-^)f{Vg{yo) + Hess5(yo)(y" - Vof + 0(' 
uj'iy'^ -yo)Vg{yofVg{yo) + 0{l) 

where the nonlinearity f2 + 0{l) again only manifests as a slow force, which is asymptot- 
ically dominated by the linear term that leads to periodic oscillations. Analogous to the 
above linear case, y^ is hence locally ergodic, the Lagrange multiplier A is asymptotically 
well-defined, and the solution x'^, y^ F-converges to the effective solution x, y. □ 

With these two lemmas, we can prove Theorem ll.2l Figured] visualizes the notations 
used in the proof to help understand the geometry. 

Sketch of the proof of Theorem Since g{q^) is at most 0{l/uj) (Lemma 13. ip . is 
close to the constrained manifold (^^^(0) in the sense that if we define for all t: 

qo{t):= mm \\q-q'^{t)\\ (40) 
ggg-i(O) 

then q^{t) — qo{t) = 0{1/ljj). Indeed, given that Vg has the maximum rank, Vg(go(i)) 
spans the normal section (i.e., the subspace perpendicular to the tangent subspace) of 



y'^ = f2{x^,yn- 
= /2(x-,y-)- 

= /2(x-,y-)- 
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Figure 1: The multiscale geometry of the constrained dynamics - slow and fast variables are x and y. 



the constrained manifold, in which q^{t) — qo{t) also lies. Moreover, g restricted to each 
normal section is an isomorphism, and both the restricted map and its inverse have 
bounded norms due to the boundedness of q'^ (i.e., compactness of the solution space) 
— this is why = 0{l/uj) implies g'^(t) - qo{t) = 0{l/u)). 

The idea is that since q^ is close enough, the constrained manifold can be locally 
viewed as a flat subspace, and F-convergence for this case has already been proved in 
Lemma 13.21 Mathematically speaking, there exists a linear isomorphism such 
that 

^,oW(9"(*)-9o(t))= [°] (41) 

where y is a vector with codimension of g^^{0) and is a null vector with dimension of 
9-^0). 

For q^{s) with s — t = 0{l/uj), we will have a full-dimensional representation: 



^,oW(9"(s)-'?o(t)) 



(42) 



and X and y will respectively be the slow and fast variables, representing the constrained 
dynamics and fluctuations away from the constrained manifold (analogous to Lemma 
13.21). This is because 



A^^^^,^{-VV{q-{s))-u:'g{q-{s))Vg{q-{s))) 



fi{x,y) 



+ 



0(1) 

-uj^giy)Vg{y) + 0{l) 



(43) 



where fi and /2 are defined as Aqy(j)(— Vy((7'^(s)). The top 0(1) on the right hand side 
of (|13]) is because ^go(j) rotates the normal section to the y-direction, i.e., 



^,o(t)V5('Z"(5)) = ^,oW(V5(go) + 0{1M) 



+ 0{l/uj) 



(44) 
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where * is some non-zero expression, and certainly 0{1/lj) = 0(1). 

The bottom on the right hand side of (03]) can also be intuitively obtained by using 
an analogous geometric argument, together with the Taylor expansion. 

Since ()^3]) corresponds to the locally flat system ([33]) . Lemma ()3.2p proved the ex- 
istence of an equivalent Lagrange multiplier as well as the F-convergence towards it. 
Moreover, ()43p and the global dynamics near the curved constrained manifold ([7]) is 
linked via a coordinate transformation i— )• Aq„((7'^ — qq), which, naturally, is slowly 
varying as qq changes. Since averaging via F-convergence (Theorem 1.2 in [31]) still 
works if the slow and fast variables are images of the original variable under a slowly 
varying diffeomorphism, the global dynamics ([7]) is F-convergent to a solution of Q. 
Notice that g{q{t)) = in ([9]) is automatically satisfied, because \m\^^^oo g{(f {i)) = 0. 

Lastly, the solution to ^ is also the solution to This is by the existence and 
uniqueness of the solution to differential algebraic equations with initial conditions. 
■ 

Remark 3.1. The above proofs show that g{q^{s)) is not only 0(tj~^) but in fact 
0(a;-2). 

Remark 3.2. li — MvHt^^ao^'^ g{(t^ {t)) exists, then ([8]) simplifies to 

\{t) = - hm u^giq'^it)), (45) 

and the solution of the modified system ([7]) will not only F-converge but also strongly 
converge to that of ([9]). Such a strong equivalent Lagrange multiplier, however, does not 
always exist, at least it is not observed in our numerical experiments nor in analytical 
investigations of linear systems (results not shown) . 

An equivalent Lagrange multiplier in the sense of averaging ([H]), as shown above, 
always exists. For numerical demonstrations, see Figure [3] (double pendulum) or Figure 
[TU] (water cluster). 

Remark 3.3. We only provided the main lines of the proof of Theorem 11.21 The details 
pose no difficulty, and are analogous to the analysis done in [311 [32] . We chose not to 
do so because full details would add 20 more pages to this paper and we do not want 
to diffuse the attention of the reader from the key contribution of this paper: SyLiPN's 
symplecticity (although the authors have not seen a mathematical justification of the 
modified mechanical system elsewhere, the modified system per se is, again, already 
widely acknowledged). 

4 Examples 

4.1 Double pendulum 

Implementation: Consider a double pendulum system. One way to represent the 
system is to use 4 degrees of freedom and 2 nonlinear constraints in Euclidian coordinates. 
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Written in the above notations, we have 



M 

V{xi,yi,X2,y2) 
9{xi,yi,X2,y2) 



mi 














mi 














m2 














m2 



-gyi - gy2 



x\ + yl 



{x2 - xiY + (y2 - yiY 



^2. 



(46) 

(47) 
(48) 



where mi,m2 are two masses, g is the gravitational constant, and Li,L2 are lengths of 
the two pendulums. To simplify our notations, we adopt a dimensionless convention and 
assume mi = m2 = g = 1. 

Implementation of our constraint-free approach on this system is straightforward 
(Eq. m). 

To integrate in generalized coordinates 6,(1), let xi = Lism6,yi = —Licos9,X2 = 
Ll sin 9 + L2 sin cj), 2/2 = —Li cos 9 — L2 cos (j). Then the Lagrangian L{q, q) = ^qMq^ — 
V{q) with q = [xi yi X2 ^2] simplifies to be 



L{9,4>,9,4>) 



1 



+ 2LiL2(cos 9 cos ((> + sin 9 sin 



+ 2Li cos 9 + L2 cos 4> 
(49) 

in which the length constraints are intrinsically handled. Corresponding Euler-Lagrangian 
equations will give the constrained dynamics. Numerically, one approximates the ac- 
tion J^'^'^^^^ L{9,(j),9,<j))dt using a quadrature rule and obtains a discrete Lagrangian 
Ld{9k,4'k,(^k+i,4'k+i)- Applying the least action principle again, a set of discrete Euler- 
Lagrangian equations [22j are obtained. Notice that the mass matrix in the generalized 
coordinates is 



Mi9, 



2Ll 



LiL2{cos 9 cos <j) + sin sin ( 



LiL2(cos ( 



cos I 



• + sin 9 sin ( 



7^2 

^2 



(50) 



which is no longer constant but position dependent. As a consequence, variational inte- 
grators given by discrete Euler-Lagrangian equations, even symplectic Euler or Verlet, 
will be implicit. 

Regarding numerical approximations to the continuous Lagrange multiplier system 
([6]), the well known algorithms of SHAKE and RATTLE can be viewed as variational 
integrators with Ist-order and 2nd-order quadrature discretizations of the Lagrange- 
d'Alembert principle f (-L(xi, yi, X2, 2/2) + ^ - gixi,yi, X2,y2))dt, in which constraints are 
realized via a vectorial Lagrange multiplier A(t) [22]. These methods also involve implicit 
solve at each step to compute the virtual force (A) that reinforces the constraints. 

Importantly, this system is much less trivial than it might appear, and the reason 
lies in the nonlinearity of the constraint g. In the modified system, for example, the 
additional constraining potential uj'^g^ g/2 is not quadratic, and it results in cubically 
nonlinear ODEs. 
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Figure 2: Integration errors by big step implicit SHAKE, small step explicit Variational Euler on 
the modified system (Eq. [7|, big step implicit Newmark (with nonlinear systems fully solved at each 
timestep), big step linearly implicit Newmark (with the nonlinear system at each step solved to the 
first iteration — no longer symplectic), and big step SyLiPN. The benchmark is provided by small 
step Variational Euler in generalized coordinates (implicit). Initial positions are 11(0) = 0,i/i(0) = 
— 1, 0:2(0) = 1, 2/2(0) — —2 and initial momenta are zero, L\ — 1 and L2 = \/2, and total simulation time 
is 100. The modified system (corresponding to Rows 3-6) uses uj = 20. SyLiPN uses (3 — 0.4. Row 3 
uses h = 0.1 /lj for stability, and Row 2,4,5,6 use a lOx bigger h — 0.05. Position errors are only shown 
on X2 and j/2 for readability. 
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Modified mechanical system 

0| < < < 




Time 

Figure 3: Comparison of Lagrange multipliers computed by ((S)) from the modified mechanical system 
(simulated by SyLiPN) and by SIfAKE. Here Lagrange multiplier at each time is a 2-dim vector (in- 
dicated respectively by red and blue). All experimental parameters are the same as in Figure [J] The 
integral in (|SJ is numerically approximated by an empirical average at time points t, t + O.l/o;, ■ • • , 
t-|-0.2-0.1/cj, t + 0.2. 

Results: Figure [2] provides an illustration of the numerical errors of different integra- 
tion methods. Symplectic numerical methods based on the constraint-free modified sys- 
tem (Row 3,4,6) produce small errors comparable to that by standard iterative SHAKE 
(Row 2). 

The original Newmark (Integrator [2]) with only the first iteration taken when solving 
the nonlinear system (Row 5) also integrates the modified system, but its error is very 
big. This is because executing only the first Newton iteration is equivalent to linearizing 
the nonlinear update equations (analogous to how we obtain SyLiPN from Push-forward 
Newmark), which can be shown to be not symplectic any more. Since it can also be 
shown that linearized Newmark is a convergent method (i.e., with finite time simulation 
error arbitrarily small by choosing h small enough), its unsatisfactory performance is 
completely due to the lack of symplecticity — notice that the energy drift is drastic 
because the system is stiff. In the nonlinear solving perspective, this is equivalent to 
saying that the predictor from the last step does not give an accurate enough initial 
guess, and Newton solver has to take more than one iteration to solve the updating 
equation in order to keep the integration symplectic. 

On the contrary, SyLiPN (Row 6) yields small numerical errors that are almost 
identical to those by the original Newmark (with the nonlinear system fully solved at 
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each timestep; Row 4). This is due to the nontrivial choice of Push-forward Newmark as 
the method to linearize, whose symplecticity survives the hnearization. At the same time, 
it is not surprising that SyLiPN produces errors shghtly larger than that by SHAKE, 
because when to is finite, the modified system is just an approximation to the constrained 
dynamics, and this approximation induces additional error besides the integration error; 
Row 3 also indicates such an approximation error, and an additional discussion on it can 
be found in Section W?2^ . 

Figure [3] compares the Lagrange multiplier, computed by SHAKE with the equivalent 
Lagrange multiplier obtained from the modified mechanical system ([7]) via averaging 
([5]). The agreement between their trajectories validates numerically the equivalency 
between the well-established Lagrange multiplier approach ^ and the modified system 
([7]) (Theorem ll.2t notice that the constraint g is not an affine function but nonlinear). 

Speed wise, generalized coordinate implicit VE (benchmark), SHAKE, Variational 
Euler, Newmark with full nonlinear solve, Newmark with one-step nonlinear solve, and 
SyLiPN respectively spent 77.9, 5.8, 0.8, 8.1, 0.3 (4.6 if HessF is not analytically provided 
but approximated by the nonlinear solver), and 0.3 seconds on the above simulation (on 
a 2.4GHz laptop running MATLAB 7.7 and 'fsolve' as the nonlinear solver with its 
default 10~^ termination tolerance). Except VE on the modified system, all methods 
use a large timestep independent of uj. Generalized coordinate approach, SHAKE, and 
Newmark are based on solving nonlinear equations, which significantly slowed down the 
computation. Linearized implicit integrators use a large step and only involve solving 
linear systems, and therefore are superior in terms of computational efficiency. 

Numerical verification of g(q)=0(co"^) 

2| , , , , 1 




20 40 60 80 100 

(0 



Figure 4: aj^||g(g(T))|| numerically computed by SyLiPN as a function of u). u) is enumerated from 10 
to 100 with an increment of 0.1, T = 50 is fixed, and experimental parameters are the same as in Figure 
[21 except that a smaller timestep h — 0.01 is used for a more accurate demonstration (the timestep is 
still coarse, because Variational Euler or Velocity Verlet requires h — 0.001 for stability). 

In addition, Figure H] numerically shows that a;^||(7(g(T))|| = 0(1) by co values, and 
therefore verifies that g{q) = 0{l/uj'^) (Remark 13. ip . 

Also worth mentioning is, given the same integration step size, a linearly-implicit 
method (e.g., SyLiPN) is not necessarily slower than an explicit one with complicated 
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force evaluations. In fact, we analytically pre-compute a position dependent nonlinear 
function, which corresponds to the inversion of the matrix I + Hessy(x)/3/i^, so that 
SyLiPN (Integrator [1]) becomes entirely explicit. Doing so, however, does not result in 
a gain in speed, because it is expensive to evaluate this pre-computed matrix function, 
which is then multiplied by a vector to accomplish the force evaluation. In fact, we 
counted the time elapses of one implementation based on linear solves and another based 
on explicit multiplications, and they respectively spent 0.3 and 4.2 seconds (the inverse 
matrix was generated automatically by Mathematica function 'Inverse', and therefore 
the code for computing it may not be optimized). In this sense, a linearly-implicit 
method is good enough. 

4.2 Convergence test 

By Corollary 12.11 and Remark 12.61 we know that SyLiPN is a 2nd order integrator. How- 
ever, this is only a quantification of its integration ability at small h. In this subsection, 
we will numerically investigate two additional questions: 

(i) How does the integration error of SyLiPN scale with the timestep h when h is 
large? 

(ii) Besides the error of ODE integration, the constrained dynamics model ([7]) is also 
just an approximation. Combining both approximations, how does the error of SyLiPN 
simulation of the constrained dynamics depend on h? 

Double pendulum: First, continue our numerical experiments on the double pendu- 
lum system (Section 14. ip . The same experimental setting will be used, but we decrease 
the simulation time from 50 to 10 so that the following investigation can be done within 
a short time. In addition, we further decrease h from 0.005 to 0.0001 in the generalized 
coordinate implicit VE so that the benchmark is more accurate. 

To investigate the convergence rate, we compute log^, error (a/i)/error(/i) as a function 
of h, where error(/i) is the numerical integration error when the integration step is h, 
and a is a constant close to 1 (we chose a = 0.9). This should converge to 2 for a 
second-order method. More generally, assuming 

error(/i) = Ch^ (51) 

for some constants C > and p, it immediately follows that 

p = log„ error(Q/i)/error(/i) (52) 

To study question (i), we obtain error (/i) by comparing the SyLiPN integration of ([7|) 
with timestep h and a benchmark integration of the same equation obtained by Velocity 
Verlet with a very small timestep COOOl/w. As can be seen from Figure El the conver- 
gence rate starts with being faster than cubic when h is macroscopic (independent of 
l/o;), but eventually approaches being quadratic when h is microscopic (i.e., h <C l/w, 
1/uj = 0.05 in our case), which is consistent with Corollary 12.11 Moreover, peculiar 
behavior exhibits in the mesoscopic region (i.e., h ~ 1/w), where the convergence rate 



19 



Rate of convergence 



03 



O 
CO 

< 



0.00000032 



Rate of convergence (zoomed-in y-axis) 
4 




Figure 5: Convergence rate and integration error of SyLiPN as functions of h. h is 
enumerated as 0.0001, 0.0002, • • • , 0.0099, 0.01, 0.011, • • • , 0.099, 0.1, 0.11, • • • , 0.2. 



frequently changes back and forth between large positives and small negatives. This 
phenomenon, together with the super-quadratic fast convergence at macroscopic h, ex- 
plains why the error does not blow up when h ^ 1/uj and agrees with the stability of 
SyLiPN. 

To investigate question (ii), we carry out a similar analysis, but this time the bench- 
mark is obtained by generalized coordinate implicit Variational Euler with a very small 
step h = 0.0001. As can be seen in Figure [H the convergence rate this time decreases 
from super-quadratic (at macroscopic h) to (i.e., non-converging, at microscopic h). 
This is because the modified system d?]) introduces an 0{l/uj) error from the true con- 
strained dynamics ([5]). This error could be reduced by increasing u (plot not shown), 
but a very large u (and a corresponding small h; also notice that the true parameter 
in the equation is w^) is rarely necessary, unless a high precision beyond 1/100 is in 
demand. 

In fact, similar /i-dependence of the convergence rate is observed when uj = 100, but 
we will not repeat the results. We chose to show plots for w = 20 because this is the 
previously used value (in Section 14. ip . 



20 



Rate of convergence 




Figure 6: Convergence rate and error of SyLiPN when compared to the constrained 
dynamics, h is enumerated as 0.001, 0.0011, • • • , 0.0099, 0.01, 0.011, • • • , 0.099, 0.1, 
0.11, 0.2. 

Moreover, we see in Figure [7| that the numerical error of SHAKE actually grows 
back bigger and then stabilizes as h is further decreased to very small values. At the 
same time, SyLiPN is not convergent to the benchmark either, but its error is one- 
order-of-magnitude smaller. This non-vanishing error of SHAKE is due to the inevitable 
numerical errors in solving nonlinear systems: we plotted the Lagrange multiplier A in 
SHAKE as a function of time (results not shown) and found that A was in the first 
few steps but then suddenly became very large; this is because h is so small that in the 
first few SHAKE steps the nonlinear solver thinks the constraints are satisfied without 
the need for any virtue force, but this leads to that the constraints are more and more 
poorly satisfied, and eventually a large Lagrange multiplier suddenly has to be chosen, 
which results in a non-zero numerical integration error. 

One might tend to attribute this issue of SHAKE to an inappropriate choice of the 
termination tolerance used by the nonlinear solver; this, however, is not true. In fact, 
no matter how small the termination tolerance is, there always exist a small enough ho, 
such that the numerical error of SHAKE will not be further reduced by choosing h < ho. 
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Error ^ Error 




(a) SHAKE (b) SyLiPN with modified system 



Figure 7: Errors of SHAKE and SyLiPN simulations of tlie constrained dynamics as functions of h. 
h is chosen from 10~^ to 10""^ with an increment of 10~^. There is no uj involved in SHAKE, but for 
comparison with SyiLiPN we mention that 10~^ <^ ho/i-o and 10~^ < ho/u. 

For instance, Figure [7] was obtained by using 'fsolve' in MATLAB as the nonlinear solver 
with termination tolerances on function value and variable both 10~^ (its default); if we 
decrease these tolerances to smaller values, for instance 10~^, smaller /I's (e.g., 10"'' as 
we tried) will make the final time integration error big. Some other may suspect the 
inaccuracy of the benchmark, but this is not the reason either (see below). 

Uniform circular motion: To rule out the possibility of an inaccurate benchmark, we 
repeated the above experiments on a particle with no potential energy but constrained 
by + = 1, whose dynamics can be solved exactly and corresponds to a uniform 
circular motion. We obtained similar results (with slightly different numerical values 
of powers), and SHAKE still exhibits inaccuracy at tiny /i, but we will not repeat the 
details. 

4.3 High dimensional case: a chain of many pendulums 

Now consider a high dimensional generalization: a chain of finitely many pendulums 
(which approximates a rope, except ropes in reality are subject to dissipations and 
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therefore not chaotic). The system is similarly modeled by Eq. [7] with: 



M 



V{xi,yi, ■ ■ ■ ,Xn,yn) 



5 yi 1 ■ ■ ■ ) ^nj Un) 



mi 
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m„ 
rUn 



i=l 



xl + yl- Li 



{x2 - xif + {y2 - yif - Ll 



{Xn - Xn-lf + (yn " Vn-lY " ^4. 



(53) 



(54) 



(55) 



where n indicates the total number of pendulums. We again assume without loss of 
generality that rrii = 1 and g = 1. 



SHAKE h=0.02 




Figure 8: Comparison of xi(t),yi(t),Xn{t),y„{t) integrated by implicit SHAKE with h = 0.005, h = 
0.01 and h — 0.02 and lineariy-implicit SyLiPN with h = 0.01 on the modified system (Eq. [T]). n — 20; 
initial positions are Xi{0) = i, yi{0) = for i = 1,2, ... ,n and initial momenta are zero; Li — 1; the 
modified system uses uj = 100; total simulation time is 20; SyLiPN uses /3 = 0.4. 



Figure [8] provides a comparison between SHAKE and SyLiPN. The system is chaotic, 
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meaning that even the same converging integrator with shghtly different timestep lengths 
win eventuahy produce completely different trajectories, and therefore we terminate the 
integration before chaotic behavior starts to manifest so that the comparison still makes 
sense. Such a termination time is decided so that SHAKE with different integration 
step lengths produce the same trajectory, but simulations beyond that time will yield 
significant deviations. SyLiPN agrees well with SHAKE till this termination time. Notice 
that we only present trajectories in Figure [8] but not integration errors. This is because a 
benchmark cannot be provided due to the chaotic nature of the system and the lack of an 
analytical solution; furthermore, the generalized coordinate approach is not applicable 
any more, because the mass matrix (high dimensional version of (|50p ) will be dense and 
complicated, and corresponding nonlinear updating equations will be very difficult to 
solve. 

Speed wise, SHAKE with h = 0.01, h = 0.005, h = 0.02 and SyLiPN with h = 0.01 
respectively spent 18.9, 37.1, 11.5 and 1.5 seconds on the above simulation (on a 2.4GHz 
laptop running MATLAB 7.7 and 'fsolve' as the nonlinear solver with its default 10~^ 
termination tolerance). Again, SyLiPN based on linear solves demonstrates a clear speed 
advantage. 

4.4 Molecular dynamics of water cluster 

We also use our method to simulate molecular dynamics. The demonstration example 
here is a cluster of water molecules, each interacting with others via non-bonded inter- 
actions of electrostatic and van der Waals forces, both of which are highly nonlinear. 

Water model: The detailed potential that we use is the one in the prevailing TIP3P 
model (see for instance |14j)). A system with N molecules has 9N degrees of freedom, 
because each water molecule consists of one oxygen atom and two hydrogen atoms, 
and each atom in 3D space is represented by 3 degrees of freedom. Indicating the ath 
molecule's ith. atom by position qai and momentum pai (each is three-dimensional), we 
write the Hamiltonian of the system as: 

« = E E ipS^r v.. + E E (e E ^ + ^n^- - ^i^] i^) 

a=li=l a=l b=a+l \i=l j=l """"^^ «2,62 ^2,62^ 

where rai^bj '■= Wlai ~Qbj\\ is the inter-atom distance, mi = correspond to the mass of 
the hydrogen atom and m2 the oxygen atom, Kc is the electrostatic constant, Qi is the 
partial charge of atom i relative to the charge of the electron, A and C are Lennard- Jones 
constants that together approximate van der Waals forces. 

In this TIP3P model (as well as most other prevailing models such as SPC, BE, 
TIPS2, and TIP4P, which are discussed in, e.g., [2], [M], or references therein), each 
water molecule is considered as a rigid body, with lengths of the two 0-H bonds and the 
H-O-H bond angle fixed as constants roH and chhoh- Detailed values of these model 
parameters could be found in, for instance, [H] . 
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Therefore, the following vectorial constraint (3A^-dimensional) will enforce the ge- 
ometry of the water molecule: 



{qu - qi2){qu - qi2f 
{qi3 - 912) (gi3 - 912)^ 

{qn - 913) (gii - ^is)"^ 



- r, 




9{q) 



(57) 
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where vhh '■= '^'i~OH sin^anoH f^) is again a constant. Consequently, the modified system 
(which SyLiPN simulates) is modeled by the following Hamiltonian: 



n = n + ^u^^Yl ((^'1,-2 - rhn? + {rl3,a2 - rlu? + (^ai,a3 " rluf) (58) 



Notice that the constrained problem is highly nontrivial (even though we already recasted 
the bond angel constraint as an equivalent constraint on the H-H distance), because the 
penalty terms are quartic nonlinear functions but not quadratic. 

Constant temperature simulations: A water cluster at a constant temperature is 
the most frequently simulated case, because (1) thermal fluctuations are indispensable 
components when modeling the reality, and (2) the N-body system is chaotic in nature, 
and therefore its long-time deterministic simulation has limited applicability. Based on 
a purely demonstrational purpose (this might not be the best model), we adopt a well- 
established approach of Langevin (see for instance [28]) to simulate the water cluster 
in a constant temperature environment. In this model, the molecules are perturbed by 
environmental noise, and such perturbations gradually balance with the internal dissipa- 
tions so that thermal equilibrium is eventually achieved. Mathematically, we integrate 
the following Langevin stochastic differential equations that model the stochastic version 
of the modified mechanical system: 



where is a 9A^-dimensional Brownian motion, (3 is the constant temperature, and 
7 is a positive constant indicating the strength of the dissipation. The system admits 
an invariant measure, which is the constant temperature Boltzmann-Gibbs distribution 
(BG for short; also known as the canonical ensemble), given by: 



N 



a=l 




(59) 




(60) 



where Z = 



fl8N 
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To simulate (|59|). we use the approach of Geometric Langevin Algorithm (GLA for 
short; see [1]), for GLA is shown to produce trajectories that are not only path-wise 
accurate but also convergent to BG ([60]) under reasonable assumptions. Specifically, we 
compose the one-step update of SyLiPN ([T|) with the exact flow map of an Ornstein- 
Uhlenbeck process (OU for short, given by dp = --fM-^pdt + ^/^q^dW) to obtain 
a stochastic version of SyLiPN. One link to the classical methods is that if a first-order 
approximation of the OU flow map is used instead of the exact flow, the composition 
will result in the well-known Euler-Maruyama scheme. 

Similarly, we compose the SHAKE one-step update with the OU flow map to extend 
SHAKE to the constant temperature setting. Under some technical conditions, the 
resulting method will approximately sample a constrained BG distribution (see Remark 
2.1 in i]): 

7r(g,p) = i-iexp(-/3i/), (61) 

where Z = Jrp,g_i^Q^exp{—l3T-l)dqdp is a different partition function, since the phase 

space is now the cotangent bundle T*g~^{0) with its base being the constrained manifold. 

One may ask about the relationship between (|60p and ()6ip . Numerical experiments 
suggest that ([5U|) converges to ([^T|) as a; — )■ oo (see the below experiment). Analytically, 
we know the equivalency between each individual constant--ff-orbits (see Theorem 1 1.2p . 
However, a full answer to this interesting yet difficult question (the proof of Theorem 
11.21 is already highly-nontrivial, let alone its thermodynamics extension) is beyond the 
scope of this paper, which focuses on the symplecticity of SyLiPN. 

On the choice of modeh As mentioned above, there are many successful water mod- 
els besides TIP3P. In addition, there are many other popular settings in which water 
molecules are coupled with the environment in difi^erent ways, including the constant 
temperature and constant pressure scenario (considered, for instance, in [3]). However, 
since the aim of this paper is to solve general problems (not limited to molecular dynam- 
ics), this section is more for an integrator-demonstration purpose rather than to claim 
any expertise in molecular dynamics; this purpose, we hope, will be well served as long 
as SyLiPN and classical constrained dynamics integrators produce agreeable results on 
nontrivial models, and therefore we did not focus on a best choice of model. 

Numerical results: In water simulations, one is often interested in the statistical dis- 
tribution of interatomic oxygen-oxygen distances in the thermal equilibrium limit, also 
known as the 00 radial distribution [14] . Figure [9] shows the histograms that approx- 
imate such a distribution, respectively obtained by long time simulations of SHAKE 
and SyLiPN. Here we chose a system with appropriately many water molecules (A^ = 7, 
63 degrees of freedom), so that peaks in the distribution could be clearly discerned. 
Since both the locations and sizes of all the peaks agree very well between SHAKE and 
SyLiPN, this is a strong numerical evidence that SyLiPN works for this complicated 
constrained stochastic system. In order to get this result, SHAKE spent 13471.75 sec- 
onds, including 12283.50 seconds on solving nonlinear systems (with 10~^ termination 



26 



SHAKE histogram 



SyLiPN histogram 



500 




500 
400 
300 
200 
100 




Li 



Figure 9: Empirical OO radial distribution in 7- water cluster obtained by long time (T = 10000) 
simulations of SHAKE and SyLiPN. 



tolerance on variable and 10~ termination tolerance on function value — the latter 
smaller because the constraint is high-dimensional and has a big Jacobian), whereas 
SyLiPN spent 1549.12 seconds, including 66.86 seconds on solving linear systems. ~9x 
speed-up is obtained by SyLiPN. 

Regarding experimental details, the initial configuration qq is obtained by minimizing 
the potential energy function V{-) starting with a random configuration (with V > 
1000) via BFGS quasi-Newton's method (see for instance [25]). However, we manually 
terminated optimization much earlier before it converges to a local minimum (V^(go) ~ 
—24 and the local minimum corresponds to ^ ~ ~57) — this way, we make sure that 
the convergence towards BG is totally due to the (correct) performance of SyLiPN and 
SHAKE but not an initial condition with high BG density. This qo is available at 
http://www.cds.caltech.edu/~intao/Papers/Water7_qIC.txt, The initial momenta 
are all zero. Both SHAKE and SyLiPN simulate one trajectory of the system till time 
10000. Judging from when the energy stops decreasing and instead starts to oscillate, the 
system equilibrates around time ~ 1000 (results not shown). We obtain the histograms 
of oxygen-oxygen distances by q's between time 5000 and 10000, and since this is much 
later than the mixing time, accuracy in approximating the distribution will be guaranteed 
due to ergodicity and law of large numbers. Histograms are drawn using 5000 bins so 
that the distribution profiles are smooth enough. The modified system uses uj = 20. 
Both SHAKE and SyLiPN use an integration step h = 0.05 (notice that GLA based 
on Velocity Verlet requires a step lOx smaller for a stable integration of the modified 
system). Dissipation coefficient 7 = 0.01, and the temperature inverse is /3 = 50. 

We also use this example to provide another numerical demonstration of Theorem 1 1.2 1 
— Figure fTOl compares the Lagrange multiplier computed by SHAKE with the equivalent 
Lagrange multiplier obtained from the modified mechanical system ([7]) via averaging 
([8]). These multipliers again agree very well (we turned off the noise and friction here; 
otherwise stochasticity will lead to diff'erent results and void the comparison). 

Although the above demonstrations already involve 63 degrees of freedom, we now 
further increase to 100 (900 degrees of freedom) to demonstrate the method's ap- 



27 



7-water-cluster: first oxygen X — modified meciianical system 




Time 

Figure 10: Comparison of Lagrange multipliers computed by ^ from the modified mechanical system 
(simulated by SyLiPN) and by SHAKE on water example. Only the Lagrange multiplier corresponding 
to the first oxygen atom is shown (3-dimensional vector). The integral in (O is numerically approximated 
by an empirical average at time points t, t + 0.1 /uj, ■ ■ ■ , t + 0.2 — 0.1 /uj, t + 0.2. Due to the complicated 
dynamics of A, lj is decreased to 500 so that visually no intrinsic details will be wiped out by the empirical 
averaging (notice that after all the equivalency is only an cj — >■ oo asymptotic result). The demonstration 
is not over long time because numerical errors accumulate rapidly in deterministic simulations of a chaotic 
system. 

plicability to large-scale problems. Figure [11] shows histograms that approximate the 
00 radial distribution in such a large system. The location and size of the first peak 
(roo ~ 2.7) respectively obtained by SHAKE and SyLiPN again agree well. Subsequent 
peaks, which are abundant, are however not easy to compare, because the simulation 
time is not long enough (although it is past the equilibration time) and hence there are 
not enough sample points to generate the empirical distributions with a small enough 
variance. Nevertheless, it is rather clear that SyLiPN produces results that agree with 
SHAKE. 

In this experiment, SHAKE spent 42234.41 seconds, including 14272.05 seconds on 
solving nonlinear systems, and ~28000 seconds on force evaluations (the computations 
of V and W), whereas SyLiPN spent 30059.47 seconds, including 318.98 seconds on 
solving linear systems, and ~29000 seconds on force evaluations (the computations of V, 
W and Hess V). Notice that, although SyLiPN still saves a lot due to the linearization, 
most of the computation has been spent on force evaluations, and therefore SyLiPN does 
not improve the simulation efficiency as significantly as before. However, many brilliant 
approaches have been proposed to greatly accelerate force evaluations for molecular 
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Figure 11: Empirical OO radial distribution in 100- water cluster obtained by simulations of SHAKE 
and SyLiPN till T = 1000. 



dynamics, including the fast multipole method [8], or simply the idea of adopting a 
cut-off distance for ignoring weak long-range forces. Nevertheless, since the topic of 
force evaluations is parallel to designing integrators, we did not employ any of these 
sophisticated methods in our simulations, but only comment that they can greatly reduce 
the cost of force evaluation, hence render the acceleration by SyLiPN more significant 
in large systems. 

It is important to notice that the times spent on force evaluations by SHAKE and 
SyLiPN are comparable; in other words, the computation of the Hessian in SyLiPN does 
not incur much computational overhead. This is because the potential is a function of 
relative distances, schematically r^j = \\xi — Xj\\. Notice that for any / 

d"^ f{r) dr d"^ f dr ^ df d'^r 
dxidxj dxi dr'^ dxj dr dxidxj ' 

but ^ and ^ are already computed while obtaining the gradient, and and q^q^ 
are not more expensive to compute. Furthermore, if a gradient method (such as Newton) 
is used as the nonlinear solver in Lagrange multiplier or generalized coordinate methods, 
the Hessian of the potential is either required or numerically estimated by the nonlinear 
solver, because the equation itself involves the gradient of the potential. In this sense, 
the computation of the Hessian is not even an overhead, unless a nonlinear solver that 
requires no gradient information is employed. 

Not only that, the linear system associated with the Hessian can also be solved effi- 
ciently. For water clusters, the Hessian is in fact strongly dominated by 9-by-9 diagonal 
blocks (there are A'^ of them). This is because the stiff penalty terms in (|58p are localized, 
i.e., although there are (9A^)(9A-|-l)/2 possible pairwise interactions, only 0{N) of them 
are stiff. Therefore, solving a linear system with coefficient matrix M + lA.essV{xk)j3h? 
can be done in 0{N) time. This will be the usual case for polymers, unless the number 
of bonds is not at the order of the number of atoms. Moreover, simplified SyLiPN for 
stiff systems (Integrator U where the stiffness is identified with w^) can be used 
for large-scale problems, and we actually use it in the above N = 100 water simulation 
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— after all, this simplified version is also convergent, symplectic, and stable. In sim- 
plified SyLiPN simulations of water clusters, the Hessian is completely block-diagonal, 
and the linear solve is executed molecule by molecule. Similarly, the nonlinear solve in 
SHAKE for enforcing the geometry of each water molecule was also carried out molecule 
by molecule; in other words, we have already used the optimal nonlinear solver in our 
implementation of SHAKE. 

In addition, in case that one is reluctant to code up the Hessian evaluation (which 
will not be difficult if the potential is a function of relative distances), s/he also has the 
choice of approximating it by finite difference and DOF + 1 evaluations of W. Due to 
Theorems 12.21 and 12. 3[ such an approximation will affect neither the convergence nor the 
symplecticity. 
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